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Abstract 

We introduce a single patch collocation method in order to compute solutions of 
initial value problems of partial differential equations whose spatial domains are 3- 
spheres. Besides the main ideas, we discuss issues related to our implementation and 
analyze numerical test applications. Our main interest lies in cosmological solutions 
of Einstein's field equations. Motivated by this, we also elaborate on problems of 
our approach for general tensorial evolution equations when certain symmetries are 
assumed. We restrict to U(l)- and Gowdy symmetry here. 

1 Introduction 

Numerical studies of initial value problems of partial differential equations on certain 
spatial domains have a long history both in basic research and in applied science, see for 
instance |6j and references therein. In particular, for applications in geometrical classical 
theories in physics, as for instance general relativity, Maxwell theory, but also Ricci flow 
(introduced e.g. in [HU]), there are interesting applications where the spatial domain is 
a 3-sphere. In general relativity, spatial 8 3 -topology plays a particularly important role 
for the standard model of cosmology based on the spatially homogeneous and isotropic 
Friedmann- Robertson- Walker solutions (see for instance [22 \ I40j). Beyond those simple 
and, at least for certain matter fields, well-understood models with high symmetry, 
there exist several outstanding open problems; of particular outstanding interest and 
motivation for this work here are the strong cosmic censorship and BKL conjecture in 
Gowdy vacuum solutions of Einstein's field equations for spatial § 3 - or S 1 x § 2 -topologies 
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It is not straight forward to deal with "non-trivial" spatial topologies, such as S 3 , 
numerically. Recall the well-known problems occurring at the coordinate axis in the case 
of standard cylindrical coordinates (p,(p,z) in M 3 . These coordinates degenerate at the 
"axis" given by p = 0. In typical equations, this has the consequence that derivatives 
with respect to the azimuthal angle 4> always come together with a factor 1 /p. If / is a 
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smooth function on R 3 , then a term like l/pd^f is well-behaved at the axis. However, 
when an equation with such terms is solved numerically, the "formal singularity" 1/ p can 
cause numerical instabilities. Applications with axial symmetry have a great history in 
all over science. Just to name examples of numerical studies of axisymmetric problems 
in general relativity, we list [19, 8, 34J. In this paper, we will not be interested in axial 
symmetry, however, it is instructive to keep it in mind for the following reason. Let us 
assume that we cover a dense subset of § 3 with one coordinate patch, for instance the 
Euler coordinates introduced below. Then it turns out that this leads, loosely speaking, 
to two "axes" on S 3 , each of which with similar properties as the axis for cylindrical 
coordinates on M 3 . Hence, we call this the "axis problem" also in the case of spatial 
S 3 -topology. 

Alternatively to this approach with one singular coordinate patch on the spatial 
domain, one can try a multipatch technique. The idea is to cover the spatial domain 
with several local regular coordinate patches. In general relativity, examples of imple- 
mentations of the multipatch technique are [381 lllj : of particular interest for our work 
here are implementations based on spectral methods in \?>2\ [35j I12j . In any case, the 
multipatch technique does not seem advantageous for the applications which we have 
in mind. First, its implementation is difficult, since one must find a stable and efficient 
way of guaranteeing the necessary communication between the local patches. Second, 
we are interested in cosmological solutions with symmetries, and in order to take full 
advantage of those, it is often a good idea to adapt the coordinates to the symmetries, 
even though this can mean that one has to deal with singular coordinate systems. We 
have decided to develop a single patch code based on the collocation method^ for the 
spatial discretization. It is general experience that such methods typically yield high 
accuracy [6]. Furthermore, a spectral single patch approach seems very natural from 
our geometric point of view, which we introduce in this paper. For the time discretiza- 
tion, we use the method of lines with a Runge Kutta integrator. Such a discretization 
technique for various spatial domains has a long history in computational physics, see 
for instance the references and examples in [6]. An alternative approach, which uses 
spectral discretization both in space and in time, has been reported on in [23]. However, 
to our knowledge, the case of spatial S 3 -topology has not been studied yet. 

Our aim is to find numerical solutions of systems of first order quasi-linear evolution 
PDEs, written schematically as 

d t f(t, x) + J2(A%f, t, x) • Vi)/(t, x) + B(t, x, f) = 0. (1) 

Here, the unknown / is a vector, the terms B are vector valued, and A 1 and Vj represent a 
collection of matrices and spatial derivative operators, as explained in more detail later. 
The spatial domain, represented by the abstract coordinates x, is § 3 . We assume in 
the following, without further notice, that the initial value problem for such a system 
is well-posed, i.e. that for any choice of initial data in a given regularity class, there 

1 In this paper we will often speak sloppily of spectral or pseudospectral methods when we mean the 
collocation method. 
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exists a unique solution / locally in time, which depends continuously on the initial 
data in a well-defined manner. In our applications, ignoring the "axis problem" above 
for the moment, the operators Vj can be thought of as the spatial partial derivatives, 
the quantities A 1 as symmetric matrices, and all coefficients depend smoothly on their 
arguments. In this case, the system is symmetric hyperbolic, and well-posedness follows 
|26} 128]. In particular, when the initial data is smooth, then the solution is smooth until 
it breaks down. We remark, that in principle our approach applies to other forms of 
hyperbolicity for Eq. ([T]), or even parabolic systems. 

We will often use the following equivalent geometric language. Namely, we say that 
we look for solutions of Eq. (pQ) on 1 x § 3 , where time t is interpreted as the canonical 
coordinate on the factor M of the manifold K x § 3 , and the t = const- hypersurfaces 
have S 3 -topology. In this paper, we will often be concerned with the case when / in 
Eq. ([I]) represents components of smooth tensor fields. Then we call Eq. ([I]) "tensorial 
equation" . 

This paper is organized as follows: Section [2] is devoted to the description and dis- 
cussion of our numerical technique. We show our geometric point of view which leads, in 
a natural way, to a single patch collocation discretization in space. We demonstrate that 
it allows a straight forward treatment of the "axis problem" mentioned above. Since 
we are interested, in particular, in tensorial equations, we discuss several related issues 
in Section 12.31 Namely, in order to express the tensor fields as collections of smooth 
functions we have decided to work with smooth global frames on § 3 . The formulation 
of tensorial equations in terms of smooth global frames on S 3 leads to certain problems 
in the presence of symmetries. Two particular classes of symmetries of our interest 
are discussed. For these sections, but indeed at many places of this paper, some back- 
ground in differential geometry would be helpful, which can be obtained for instance 
from [221 129] . In Section \2A\ we elaborate on our numerical infrastructure in general, 
and we discuss certain issues related to numerical stability which are present for spatial 
§ 3 -topology. In Section [3l we proceed as follows. First, we introduce the mathematical 
and physical background of our test application. Then we test the code in that setting 
and elaborate on numerical errors, stability and performance. Finally, in Section 0J we 
summarize and conclude. We comment briefly on possible problems when some of our 
special assumptions are dropped, point to open issues and list aspects for future work. 

2 Geometric ideas and numerical implementation 

2.1 Euler coordinates and the Euler map 

In the following, we consider the 3-sphere § 3 as the submanifold of R 4 given by 
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The Euler coordinates of S 3 , also appearing as Euler angle parametrization, etc. in the 
literature, is the coordinate patch given by 



X X 

x\ = cos — cos Ai, X2 = cos — sin Ai , 

X X 

x% = sin — cos A2, X4 = sin — sin A2, 



(3a) 



in terms of the coordinate functions x G [0, 7r], Ai,A2 G [0, 2tt[. We will rather use the 
coordinates {X1PI1P2} determined by 

Ai = (pi + P2 )/2, A 2 = ( Pl - p 2 )/2. (3b) 

The Euler coordinates smoothly cover the dense subset of §> 3 given, when the points \ = 
0, 7r are taken away. We expect that other choices of coordinates with similar properties 
are also appropriate. Although the motivation for choosing the Euler coordinates stems 
from Gowdy symmetry, as becomes clearer later, they are robust enough for more general 
cases. 

Certainly, the relations in Eq. ([3|) are well-defined for all X1PI1P2 G R- This is also 
true when we consider XiPi^Pi G (M- m °d 47r). Geometrically, the Euler coordinates 
given by Eq. ([3]) can thus be interpreted as a smooth map from the 3-torus 

T 3 := (R mod 4vr) x (E mod 4vr) x (R mod 4vr) 

to § 3 , which we call the Euler map 

$ : T 3 -» S 3 , 

/ s ( X Pi+ P2 X • Pi + P2 

{X,Pl,P2) ^ ( cos -cos — , cos -sin — , ^ 

. X Pi ~ P2 . X . Pi - P2 \ _ s3 

sin — cos , sin — sin G s> . 

2 2 ' 2 2 J 

The Euler map <E> is even a diffeomorphism when we restrict it to e.g. x G]0,7r[ and 
restrict the image correspondingly. But at the points x = 0, 7r, the inverse is not well- 
defined. Note, that in the whole paper we will often make the canonical identification 
of the isomorphic groups U(l), S , (R mod 2ir) and (R mod An), and henceforth not 
distinguish between them. However at this point, our definition of the map requires 
to stand on 47r-periodicity at least for the coordinate Xi but we come to the standard 
2-7r-periodicity in a moment. 

Let / be a smooth function on § 3 . In the following we consider / := /o$, where $ is 
the Euler map defined in Eq. @. Hence / is a smooth function on T 3 . For simplicity, we 
write / instead of /, and often make no difference between the "original function / on 
S 3 " and the "corresponding function / on T 3 " . If necessary, in order to avoid confusions, 
we sometimes say that "/ is a smooth function on T 3 originating in a smooth function 
on S 3 " . 

Motivated by this simple geometric relation given by the Euler map, our approach 
is the following one; the details are worked out in the subsequent sections. Our aim is 



4 



to solve partial differential equations with a spatial domain §. Let us suppose that all 
coefficient functions and unknown functions, which we want to solve for, in the equation 
are smooth functions on 8 3 . Furthermore, suppose that all derivative operators stem 
from smooth globally defined vector fields on § 3 . Since $ becomes a diffeomorphism, 
when we restrict it to a dense subset of S 3 , all these functions and vector fields correspond 
in a unique manner to smooth functions and vector fields on the corresponding subset 
of T 3 . Hence, we can solve the equation as if it its spatial domain were T 3 . However, 
because is not a diffeomorphism globally, we have introduced formal singularities to 
the equations, analogous to the singularities at the axis of cylindrical coordinates on M 3 
discussed above. One main conclusion in the following sections is that the hypothesis, 
that all these quantities on T 3 originate in smooth quantities on § 3 , allow to regularize 
the formally singular behavior at those points. By all this, we will successfully make "S 3 
periodic in all three spatial directions" , and can henceforth use spectral methods based 
on the standard Fourier basis for the spatial discretization of the equations. 

In all of what follows, we will assume, for simplicity, that all functions involved do 
not depend on the coordinate p2- We call such functions U(l)-symmetric, and later 
we interpret this symmetry geometrically. The generalization of the following ideas is, 
however, straight forward. This assumption has the following nice property. Let / be 
a smooth U(l)-symmetric function on § 3 . Then / o $ is 27r-periodic (instead of Att- 
periodic) in \ an d pi- This is so, because we can use the symmetry in p2 to switch the 
signs of all terms in Eq. @ consistently as needed. Hence, in the following, we only 
need to deal with the standard 27r-periodicity. 

Note, that analogous spectral approaches for equations with spatial domains dif- 
feomorphic to S 2 have been implemented before, see for instance (2j [6] and references 
therein. 

2.2 Analysis of Fourier series of given smooth functions on § 3 

Let / be a given smooth U(l)-symmetric function on § 3 . The corresponding function 
/ on T 3 is smooth, and hence has a representation in terms of a Fourier series. We 
come back to the question of convergence of such a series later. At this stage, we 
want to assume that this Fourier series is finite, and, without loss of generality, with N 
summands both in x an d p\- Since / is U(l)-symmetric function, it is 2-7r-periodic in 
both x and pi, as argued before. Let / be real valued. Then it must be of the form 

N 

f(x,Pi)=Fo(x) + 2Re^2F p ( X )e i ^. (5) 

P =i 
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Here, Fq(x) 1S a rea l valued function, and F p (x) is allowed to be complex for all p > 1. In 
[3J, we demonstrated that the properties of the Euler map under U(l)-symmetry implies 



( N 



F P (X) = { 



2 ^2 f' n <P cos n X + for p > even, 
n=l N (Oh) 

•2* J^/n 

iP sin nx for p > odd, 



v. n=l 

for some, in general, complex coefficients f np ; only the coefficients for p = must be 
real. Some of the factors in these expressions are chosen for later convenience. Note, 
that in [3J, the coordinates (x,pi,P2) are defined slightly differently. For example, the 
function x here must be substituted by 2% to compare to the expressions in [3]. In any 
case, we further find that for even p > 0, there are the "compatibility conditions" 

N N 

hp + 2 hn,p = 0, £ / 2 n-l,p = 0. (6b) 
n=l n=l 

These originate in the fact that the functions F p (x) must vanish at the degenerated places 
X = 0, 7T for all p > 0. For odd p, this is implied by the expression Eq. (|6a|) automatically, 
and hence there are no compatibility conditions. The corresponding expressions are more 
complicated, but analogous, when we give up U(l)-symmetry. 

Now, we want to study what happens to the Fourier series of /, when it is differen- 
tiated along a smooth tangent vector field on S 3 . By this, we mean that the abstract 
derivative operator Vj in Eq. (pQ) is of the form Vj/ = V a d x af for a smooth tangent 
vector field V = V a d x a on S 3 . Here, our convention is that x a represents three abstract 
spatial coordinates, which neither need be the Cartesian coordinates, nor the Euler co- 
ordinates used before. However, we will restrict to Euler coordinates in the following. 
Furthermore, we assume Einstein's summation convention. The coefficients V a are just 
functions on S 3 which are not necessarily smooth when we consider the Euler coordi- 
nates. Of particular importance will be the following tangent vector fields, whose origin 
we explain later and which, with respect to the Euler coordinates, take the form 

Yi = 2 sin pi d x + 2 cos p\ (cot x dpi ~ csc X d P2 ) , (7a) 
I2 = 2 cos pi d x — 2 sin p\ (cot x dpi ~ csc X d P2 ) , (7b) 
Y 3 = 28 P1 . (7c) 

The factors 2 are chosen for consistency with our discussion in [3]. Now, as we explain 
later, any smooth vector field V on S 3 can be written as a linear combination V = 
V a Y a with V , V 2 , V s smooth functions on § 3 . Hence, under the assumption that all 
differential operators in the equations stem from smooth vector fields on S 3 , it follows 
for U(l)-symmetry, that all "formally singular" differential operators in our equations 
must be of the form — F(x, Pi) c °t X^pu with F some smooth function on S 3 . Without 
the assumption of U(l)-symmetry, there can additionally be singular operators of the 
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form (F(x, PI1P2)/ si n x)d P2 ■ This shows that the formally singular terms here are of the 
same type as in the "axis problem" . The differences to the case of cylindrical coordinates 
on M 3 are twofold. First, in the case of § 3 , we have two such "axes" simultaneously at 
X = and ir. Second, the axis itself is not topologically a line here, but a closed circle. 
In the following, we restrict our attention to the operator relevant to U(l)-symmetry 
- cot xd pi ■ 

So, let / be as before, and g := d Pl f. Since g is again a smooth U(l)-symmetric 
function on §> 3 with finite Fourier series, the analogue of Eqs. §5§ and ([6]) holds with 
F p (x) substituted schematically by the function G p {x)- Then, for x 7^ 0,7r, we get, 

- cot x Gpix) 



k=l 



2i 



c 2 , p sin x + ^{ (h,p + h+l, P ) sin 2kx 

+ (c k+ i, p + c fc+ 2,p) sin(2A; + l)x} 

bl >P + 5Z{ ^ Cr 'P + Cf, + 1 -p) cos ( 2r - i)x 
+ (b r , p + br+i,p) cos 2rx| 



for p > even, 



(8a) 



N 



r=l 



for p > odd. 



Here, we define 



N 



N 



b r, P '■= 92n >P' Cr 'P := 92n-i, P for r > 1. 



(8b) 



The computations leading to this result are described in [3]. This result means that, as 
soon as the Fourier coefficients of d pi f are known, the Fourier coefficients of the complete 
"formally singular term" — cotxd Pl f can be computed. 

Now, let us consider the general case of a given smooth U(l)-symmetric function / 
on § 3 . The associated function on T 3 has an infinite Fourier representation in a general, 
with rapidly decreasing coefficients f np . This is a standard result from Fourier analysis, 
which can be found in |37} 17]. This last property means that the modules of the Fourier 
coefficients is bounded by a uniform constant times any negative integer power of the 
two summation indices n and p. This property is often referred to as "exponential 
convergence" . It is a general fact under our conditions that the Fourier series converges 
pointwise absolutely and even uniformly. We find straight forwardly that the Fourier 
series of / must be of the form given by Eqs. ©, setting N — > 00. The infinite series 
of the compatibility condition Eq. (|6bj) converges because the coefficients are rapidly 
decreasing. 

Now consider the inverse question. Let a function be given on T 3 , of the form of 
Eqs. © and ([6]), for N = 00 with rapidly decreasing coefficients f n>p . The standard 
theory implies that the series converges pointwise absolutely and uniformly to a smooth 
function / on T 3 . However, does / originate in a smooth function on S 3 ? In general 
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the answer is no, because the compatibility conditions Eq. (|6b|) are just necessary, but 
not sufficient for smoothness. Indeed, it is sufficient that for any p, the function F p (x) 
is a smooth 27r-periodic function on R, which has a zero of order p at x = and it. 
In particular, each function F p {x) e%ppl then originates in a smooth function on S 3 . The 
argument for proving that this implies that / is a smooth function on 8 3 , uses the theorem 
about Fourier series on S 3 in [37J. Namely one can show that under these assumptions, 
/ as a function on § 3 can be represented as an infinite series of spin-weighted spherical 
harmonics with again rapidly decreasing coefficients. 

Consider the derivative g = d Pl f. In particular, the formula for — cot X9 in Eqs. (JSj) 
also holds in the limit N — > oo, and the series expressions there converge at least point- 
wise at all x 7^ 0) 71 "- The function — cotxd pi f is a smooth function on T 3 because 
each F p (x) in Eq. (0) is a smooth 2-7r-periodic function in x with appropriate zeros at 
the "singular places" x = 0> 7r - This is a nice result because it shows that Eq. ([8]) is 
meaningful at the singular locations, and hence allows to evaluate the formally singu- 
lar term — cotxd pi f explicitly there, even in the limit N — > oo. However, we remark 
that — cotxdpif does not originate in a smooth function on § 3 , because Eqs. (|8|) is not 
consistent with Eq. ([6]). This is not a problem because the formally singular operator 
is only a part of a differential operator defined by a smooth vector field on S 3 . Indeed, 
the result, when this "full" differential operator is applied to a smooth function, yields 
a smooth function on § 3 . 

2.3 Symmetry and related issues for tensorial equations 

Before we discuss, how these results can be applied in practice, let us first consider some 
consequences for tensorial equations. 

2.3.1 Smooth frames on S 3 

Recall, that one of our main assumptions is that, at any given instance of time, all 
unknowns and coefficients in the equations are smooth functions, and that all differential 
operators are determined by smooth globally defined vector fields on S 3 . However, in 
order to turn tensorial equations into partial differential equations for smooth scalar 
functions, we need to introduce smooth frames on § 3 . We would like to mention that 
an alternative way of treating tensorial equations on S 3 in the case of Gowdy symmetry, 
see below, can be found in [18] for the case of spatial S 1 x S 2 -topology. 

Let us recall some well-known facts. Let § 3 be given as in Eq. ([2]). Assume the 
standard matrix representation of the Lie group SU(2) [37] . The map 

^ : S° -> SU(2), (xi, x 2 , 23,24) ^ , ■ ■ 

\x 3 + 1x4 21 - 1X2 J 

is a diffeomorphism, which can be used to transport the group structure of SU(2) to 
S 3 . Hence, both SU(2) and § 3 can be considered as identical Lie groups via the map 
$ '. Thus, from the standard SU(2) group multiplication, we can define left and right 
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translation maps, 



L,R : S 3 x S 3 -> § 3 , (u,v) \—* L u (v) := uv, (u,v) ^ R u (v) := vu, 

so that L u and R u are diffeomorphisms §> 3 —* S 3 for each point u £ S 3 . Those maps 
can be employed to construct smooth global frames. First one chooses a basis of the 
tangent space at the unit element e of the group. We choose the Pauli matrices with 
non-standard normalization 

considered as elements of T e (SU(2)). Then, one uses the push forward of L u or R u to 
transport this basis smoothly to any other point u £ SU(2) 

(Y a ) u '■= {L u )*iXa)i {Z a )u '■= (Ru)*(Y a ). 

Clearly, {Y a } is SU(2)-left invariant while {Z a } is SU(2)-right invariant and both are 
smooth global frame fields on S 3 . It is straight forward to check that they satisfy 

3 3 

[Y a ,Y b ] = 2Y t e abc Y c , [Z a ,Z b } = 2Y j e abc Z c , [Y a ,Z b }=0, (9) 

c=l c=l 

where e abc is the totally antisymmetric symbol with e 12 3 = 1. Here, the brackets denote 
the Lie bracket. For {Y a }, we have already written the explicit expressions with respect 
to the Euler coordinates in Eqs. ([7|). For {Z a }, we have 

Z\ = —2 sin p2 d x — 2 cos p2 (cot \ d pi — esc % 9 P2 ) , (10a) 
Z2 = 2 cos P2 d x — 2 sin p2 (cot \ dpi ~ csc X d P2 ) , (10b) 
Z 3 = 2d P2 . (10c) 

On Kx§ 3 with a time function t, we assume that each t = const-hypersurface is 
diffeomorphic to S 3 with Euler coordinates {x-> Pi, P2\- Hence, on each of these surfaces, 
the expressions Eqs. ([?]) and (fTUl) define the fields {Y a } and {Z a }. Geometrically, we 
thus have for all a = 1, 2, 3 

[d t ,Y a ] = [d t ,z a ] = o. 

Now we write an arbitrary globally defined smooth frame {e{\ on 1 x S 3 as follows. 
By {ej} we mean the collection of 4 frame fields {eo, e%, 62, 63}. We set 

e = d t , (11a) 

and write, 

e a = e a b Y b , (lib) 

where (e a b ) is a smooth 3 x 3-matrix valued function with non- vanishing determinant 
on S 3 . Our conventions for frames is that the index always corresponds to the "time 
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frame vector" eo and a,b,. . . = 1,2,3 correspond to the "spatial frame vectors". When 
we write indices = 0, . . . , 3, we mean both the time and spatial frame vectors. 

For a tensor field S, say for example a covariant 2-tensor, we write S%j := 5(ej,ej), 
S a b '■= S(e a ,eb) etc. We stress, that it is important to understand, that writing e^, does 
not mean the ith component of a vector e, but rather the ith vector field of the frame 

2.3.2 Symmetry reductions of tensorial equations 

Let S be an arbitrary smooth tensor field on MxS 3 . We say that S is ^-invariant, provided 
C^S = everywhere. Here denotes the Lie derivative along £. The coefficients Sij 
of 5, with respect to an arbitrary frame {e^}, are constant along £, if, and in general 
only if, [£, ej] = 0. Now, suppose that the functions Sy of such a tensor are among the 
unknowns of the system of partial differential equations which we would like to solve. 
Often, we would like to exploit the symmetry of the unknowns and reduce the equations 
to some simpler form. If the unknown functions are constant along £, such a reduction 
can be done directly, if £ has the meaning of a spatial coordinate vector field. Hence, in 
the following, we will be interested in frames such that [£, e,] = 0; in this case, we say 
that the frame is £- invariant. 

Let us consider special cases of interest for us. We have already introduced the notion 
of U(l)-symmetry for functions before. For general smooth tensor fields S on R x § 3 , we 
define it by the requirement that S is ^-invariant. One can define U(l)-symmetry more 
geometrically, but for the purpose of this paper, our definition is sufficient. The integral 
curves of Z% = d P2 are circles. Hence the symmetry group is isomorphic to U(l), which 
motivates the name. Now, it is straight forward to construct ^-invariant frames {ej}. 
Namely, for the ansatz Eqs. (|llj) . this requirement is equivalent to 

d P2 (e a b ) = 0, 

due to Eqs. ©. The consequence is, that both the frame matrix (e a b ) and the frame 
components of all tensor fields, which are ^-invariant, are constant along p2, and hence 
are U(l)-symmetric functions. Provided, our equations are formulated with respect to 
functions with this property only, the spatial domain reduces to that of the coordinates 
{x, pi}- We refer to this as the 2 + 1-reduced equations. 

Another symmetry assumption of interest is Gowdy symmetry with two spatial sym- 
metry vector fields. We say, that a tensor field is Gowdy symmetric, if it is U(l)- 
symmetric and additionally yjj-invariant. Again, this definition can be made more ge- 
ometric. In any case, let us suppose that a frame field {e^}, obeying Eqs. (fTTI) and 
[Z3,ei] = 0, is given as before. Now, it turns out, as argued in [3], that the assumption 
that {e{\ is a smooth globally defined .^-invariant frame on § 3 , does not allow that it 
is 13-invariant in addition. Hence, although we can find a frame such that the frame 
components of arbitrary Gowdy symmetric tensor fields are constant along p2, it is not 
possible to achieve that in addition they are constant along p\\ recall I3 = d pi . Thus, 
even if we assumed Gowdy symmetry, the spatial domain of our equations would not 
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reduce directly to that of the coordinate x, i.e. not to 1 + 1 dimensions. This difficulty 
is a consequence of our assumption that the frame is smooth globally on the manifold; 
if we allowed the frame to become singular at some places, then the situation would be 
different. However, we would get other problems due to the additional singularities. 

Nevertheless, even under these assumptions, it is possible to perform the following 
"indirect" reduction of the equations to 1 + 1 dimensions in the case of Gowdy symmetry. 
For the frame components of any smooth tensor field S, for example in the case of a 
covariant 2-tensor, which is Y-j-invariant, we find 

Y 3 (Sij) = %t/ + %,T/" (12) 

with 

T/ ei , := [Y 3 , ei ]. 
Under our assumptions for the frame Eqs. (jllH . we see that 

T/ = 0, Ti° = 0. 

For the spatial components, one has 

T a a ' =Y 3 (e a c )f/ + 2e a % c d f/. 

The matrix (f c a ) is defined here as the inverse of the matrix (e a c ). As soon as we fix 
the transport of the frame in time, we can compute the time derivative of the matrix 
(T a a ); we do this for a particular example in Section 13.11 In general, we can expect 
that these evolution equations for the matrix (T a a ) are non-trivial. In any case, the 
idea for the "indirect" reduction to 1 + 1 in the Gowdy case is the following: first, 
substitute the 13-derivatives of all tensor field components in the equations by means of 
Eq. (fT2|) . Second, append the evolution equations for the matrix (T a a ) to the system 
of equations and, third, evaluate the equations only at p\ = 0. Then, with respect to 
the Euler coordinates, all unknowns only depend on t and x, and the evolution system 
is closed. Note, however, that it depends on the properties of the evolution equations 
of (T a a ), whether the resulting system of evolution equations yields a well-posed initial 
value problem. In the example, which we discuss later, this is the case. 

We remark that under the assumption of U(l)-symmetry, all results obtained here 
also hold for spatial S 1 x S 2 -topology. We do not elaborate on this further; a discussion 
is given in [5]. 

2.4 Numerical implementation 

2.4.1 Discretization and our numerical infrastructure 

In order to compute approximate solutions of our system of partial differential equations 
Eq. (fT]) by means of a computer, we need to discretize the equations and the unknowns. 
Our analysis before based on Fourier series suggests spectral discretization [TJE] in space 
with the standard trigonometric basis. We follow most of the conventions in [BJ. In 
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order to keep the presentation as short as possible here, we do not write down formulas 
wherever we follow the standard conventions. In particular, we use the collocation 
method. For the spatial grid in any of the spatial dimensions, referred to as x, we set 

2-7T 

x k = (k + v) — , k = 0,...,N-l, (13) 

where N is the number of grid points in the chosen spatial direction. The quantity 
fj, E [0, 1[ is a shift quantity. For the standard collocation method which we use, the 
quantity N must be odd. 

For simplicity we use the so-called partial summation algorithm j6] for computing 
the discrete Fourier transforms (DFT) so far; however, we plan to switch to the Fast 
Fourier Transform algorithm (FFT) [101 133j . in order to optimize performance for high 
spatial resolutions. 

This discretization of the equations and unknowns in space yields a system of ordinary 
differential equation (ODE) in time for the spectral coefficients of the unknowns, or 
equivalently for the values of all unknown functions at the spatial grid points. This 
system of ODEs is called the semi-discrete system. In order to solve this numerically, 
one must discretize time as well. For this, we have implemented a couple of Runge Kutta 
(RK) variations described in [33]; namely, first, the non-adaptive 4th order RK scheme, 
second, the 4th order "double-step-adaption" RK scheme and, third, the adaptive 5th 
order "embedded" RK scheme. For those schemes, the time adaption is always "global 
in space", namely, at a given time the maximal estimated error at all spatial points is 
taken. The parameter rj is the desired accuracy, according to Eq. (16.2.7) in [33], where 
it is called Ao- The lower its value, the stronger is the tendency of the adaption scheme 
to decrease the time step h. For practical reasons, we also define a minimal time step 
hmin, so that the adaption scheme is prevented from reaching unpractically small values 
of h. 

A sophisticated discussion of errors and convergence in such discretization approaches 
is given in [7]. We will not elaborate on this, in general, very complicated problem 
analytically. Instead, we will investigate errors and convergence in our test applications 
empirically in Section 13.21 The positive experience, which people have gathered over 
many years of research with the collocation method, is summarized in Boyd's empirical 
"assumption of equal errors" [6], which we decided to rely on in our numerical work. 

It is clear that many classes of problems require adaptive techniques for the spatial 
resolution. One particular effect for underresolved numerical solutions obtained by the 
collocation method is aliasing. We have not yet implemented any of the explicit an- 
tialiasing recipes, given for instance in [6]. Instead, we use the following simple global 
spatial adaption technique so far. At each time step of the numerical evolution, the 
program computes the Fourier transform of one representative unknown function; which 
one is chosen requires some experiments. In our applications, where symmetry implies 
that only one spatial direction is significant, it is then sufficient to do the following. 
From the spectral coefficients of this unknown, the code determines the "power" of the 
upper third of the frequency spectrum with respect to the significant spatial direction, 
divided by the total power. It is straight forward to generalized this to more general 
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situations. By "power" we mean the sum over the squares of the modulus of the Fourier 
coefficients. We call the result of this computation the "adaption norm" Norm.( adapt \ 
Besides adaption itself, this "norm" can also be thought of as a measure for the aliasing 
error. When, during the numerical evolution, Norm^ adap *- ) exceeds a prescribed thresh- 
old value, the code stops automatically, interpolates all quantities to a higher spatial 
resolution and continues the run. In each of these adaption steps, we have found it to 
be reasonable to almost double the spatial resolution. In any case, note that this is a 
primitive adaption method, since it is "global in space". In particular, for solutions, 
which develop sharp localized features, a local adaption method in space would be more 
desirable. This is a future work project. However, even Gowdy spikes have been treated 
with our much simpler method in [3]. 

Let us, furthermore, mention the possibility of the Intel Fortran compiler [21] on 
Intel CPUs, which we have worked on exclusively up to now, to switch from the standard 
machine supported number representation called "double precision" with round-off errors 
of the order 10 -16 to software emulated "quad precision" with round-off errors of the 
order 10~ 32 . Indeed, this possibility is exploited in our applications, as is discussed later. 

2.4.2 Practical issues for spatial § 3 -topology 

In this paper, we discuss two codes for spatial § 3 -topology which both use the infrastruc- 
ture above, but with slight differences. For one of the codes, we assume U(l)-symmetry 
with a choice of orthonormal frame reducing the equations to two spatial dimensions; 
we call this code the 2 + 1-code. For the other code, which we call 1 + 1-code, we assume 
Gowdy symmetry, and suppose that the indirect reduction to one spatial dimension de- 
scribed before leads to a well-posed initial value formulation. The particular equations 
which we implement and study in both cases are discussed in Section 13.11 

In order to apply the results, which we obtained for the properties of Fourier series of 
smooth U(l)-symmetric functions on § 3 in Section \2. 21 in our discretization approach, we 
assume that for any choice of resolution, the solution of the corresponding discretized 
equations originates in smooth functions on S 3 . Let us suppose, that our evolution 
equations have the property, that if the initial data of the continuum problem is smooth, 
then the corresponding solution of the continuum equations is smooth. For instance, this 
is the case for symmetric hyperbolic systems. Then, if the initial data is approximated 
by functions, which originate in smooth functions on § 3 , and, if all our assumptions 
about the coefficients and derivative operators in the equations before hold, then the 
solution of the discretized originates in smooth functions on § 3 for any resolution. This 
is at least true, as far as we can neglect errors caused by aliasing and the finite number 
representation in our computer. Let us assume just for a moment, that these errors can 
be neglected. In particular, Eq. ([8]) can be applied to the 2 + 1-code directly. For the 
1 + 1-code, we need one further observation, since the equations are only evaluated at 
pi = 0, and hence there is no information about even and odd Fourier modes with respect 
to p\. Namely, due to Eq. (|6a,j) . all coefficients associated with cos-modes with respect 
to x must correspond to an even mode with respect to p\. Analogously, all coefficients 
associated with sin-modes with respect to x must correspond to an odd mode with 
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respect to p\. This information is sufficient to use Eqs. (JSJ) as in the 2 + 1-case. Now, for 
even p > 0, there are two ways of computing the coefficients b^ p and of the formally 
singular terms. Namely, due to the compatibility conditions Eqs. (f6bj) . we can write 
both 

N/2 r-1 N/2 r-1 

h N _ \ ~~ N N \ " N N _ \ N \ " N 

V r,v — / j 92n,v ~~ 2^°'P Z-/^ 2n >P' Cr >P ~~ 2—/ 92n-l,p ~ / , 92n-l,y 

n=r ~ n=l n=r n=l 

For odd p, there is only one way of writing these coefficients. Although each pair is 
equivalent in exact computations, there can be a difference numerically. We refer to the 
first way of computing these coefficients as "up-to-down" , since we need the information 
of all high frequencies to compute the low frequency coefficients recursively. The second 
variant is called "down-to-up" , since the information from low frequency coefficients is 
used to compute the high frequency coefficients recursively. 

A priori, both ways have the potential to amplify numerical instabilities. In par- 
ticular, although the solution originates in smooth functions on § 3 at one time of the 
evolution, this together with round-off errors and aliasing can cause a drift, so that the 
form given by Eqs. © are violated eventually. We have not yet built the special struc- 
ture of the Fourier series into our numerical infrastructure. Thus, it is possible, that 
such errors accumulate, such that, after some time of evolution, the numerical solution 
does not represent a smooth solution on S 3 anymore. Indeed, we found in our numerical 
experiments with the 2 + 1-code in [3], that, without precautions, the numerical solution 
typically drifts away strongly for both the up-to-down and down-to-up method. We were 
not able to pin-point the problem. However, for the down-to-up method, it turns out to 
be sufficient, after each time step, to set all those Fourier coefficients to zero explicitly, 
which are supposed to vanish according to Eq. (|6ap . With this manipulation, the numer- 
ical evolution becomes stable. In particular, Eq. (|6bp stays satisfied within reasonable 
error limits, and the code is convergent and able to reproduce explicitly known solutions. 
The up-to-down method, however, we were not able to stabilize. 

We can expect that similar practical issues exist for the 1 + 1-code. Here, however, 
we must proceed slightly differently, since the form given by Eqs. © cannot be enforced 
explicitly. In order to control the smoothness of the numerical solution nevertheless, 
the idea is to control the unknowns directly at the coordinate singularities x = 0, 7r in 
terms of "boundary conditions"!!. Let the frame {e^} be i^-invariant as before, and let 
S be one of the unknown ^-invariant tensor fields. Since I3 = 1LZ3 on the symmetry 
axes, one obtains Y%{Sij) = there. Exploiting this information by means of Eq. (|12|) . 
implies a homogeneous linear algebraic "boundary system" , which yields the "boundary 
conditions" for S 

Si'jTi + S ir T j j ' = at X = 0, tt. 

For our particular test case later, we solve the boundary system in Section 13.1.31 Let 
us assume for the moment that we have solved this system. In general, we would like 



2 This is the terminology from [18]. We use it despite the fact that there are no geometrical boundaries 
at x = 0, 7r. 
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to have the possibility of either letting the numerical evolution proceed freely and just 
monitoring, how well those boundary conditions are satisfied, or, if necessary, we would 
like to enforce the boundary conditions. The latter means that we set the values of the 
unknowns at x = 0, 7r explicitly to the values implied by the boundary conditions. In 
order to make this possible, we modify the spectral conventions slightly as follows, so 
that both boundary points x = anci X = ^ correspond to grid points. Let / be some 
unknown function, which we discretize as 



In the standard collocation approach, which we use for the 2 + 1-code, the number 
of grid points N, according to Eq. (|13p . is odd, and M = N. Recall that the dis- 
crete Fourier transform [6] is the linear map from the values of / at the grid points 
(f(x ), . . . ,f(x N -i)) to the Fourier coefficients (a Q ,bi,ai, . . . , &(m-i)/2j a {M-i)/2), which 
is bijective for these choices of N and M. This is true for any choice of shift /x. For the 
1 + 1-code now, we choose \i = 0, any even number N, and set M = N + 1. In this 
case xq = and Xjv/2 = 7r * ^ ne nna ^ s easily that for this, the standard discrete Fourier 
transform is the map 



Hence, the map has the standard properties except for the highest frequency. The 
fact, that the discrete Fourier transform always yields zero for the highest sin- mode 
can be understood easily, because the value of sinn^ is always zero for n = (M — l)/2 
at x = kty. The main point is now that this map is nevertheless invertible. For 
spectral differentiation, we ignore the frequency n = (M + l)/2 completely. In practice, 
it is expected that this is not problematic, since the highest frequencies are typically 
insignificant. Now, in our numerical experiments, which have so far restricted to the 
down-to-up method, we find that the 1 + 1-code with this spectral infrastructure is very 
stable and convergent. This is true even without enforcing the boundary conditions 
at all and hence without any explicit control of the smoothness in regimes where the 
solution is relatively smooth. Recall that this is not so for the 2 + 1-code. Furthermore, 
the violations in the boundary conditions typically converge to zero with increasing 
resolution. However, if the simulation approaches a non-smooth regime of the solution, 
it seems often necessary to enforce the boundary conditions; this is discussed for our test 
applications. 

3 Analysis of test applications 

Before we test and analyze our numerical method in Section 13.2^ we briefly introduce 
some background for our test application in Section 13.11 More details can be found in 
our similar discussion in [4], where we emphasize the physical and mathematical ideas 
and interpret the results. 




(M-l)/2 



(f(xo), ■ ■ -,f(xN-i)) >-> (ao,bi,ai, . . . ,b N/2 -i, a N / 2 -i,0;2a N/2 ). 
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3.1 Background of the test application 
3.1.1 Physical and mathematical background 

Our aim is to compute cosmological solutions of Einstein's theory of relativity; in partic- 
ular we are interested in the strong cosmic censorship conjecture in Gowdy vacuum solu- 
tions of Einstein's field equations for spatial S 3 - or S 1 x § 2 -topologies [21] 1251 l9| [18] \$\ [36] . 
All our discussions assume vacuum and a cosmological constant A, so that Einstein's field 
equations (EFE) in geometric units c = 1, G = 1/(8tt) read 

R = Xg. (14) 

Here g^ is the spacetime metric, which is the fundamental unknown encoding the infor- 
mation about the gravitational field. Its Ricci tensor R [22] is a 2nd order quasi-linear 
expression in the metric. We will always assume four spacetime dimensions, that the 
signature of the metric is Lorentzian (—,+,+,-1-), and that Cauchy surfaces, i.e. the 
"surfaces of constant time", are diffeomorphic to 8 3 . Furthermore, we suppose A > 0. 

In [30] 131] . Penrose introduced his notion of conformal compactifications. The idea 
is to rescale the physical metric g by means of a conformal factor f2, which is a smooth 
strictly positive function on the spacetime manifold M. This yields the so called confor- 
mal metric 

g := fl 2 g. 

Now, loosely speaking, if it is possible to attach those points to M, which are the limit 
points of vanishing £1, so that the new manifold M is smooth and the metric g can 
be extended as a smooth metric on M, then we say that the original spacetime has a 
smooth conformal compactification. The references above, but in particular [15], give 
further necessary technical requirements to make this loose statement rigorous. Under 
those conditions, the set f2 = is a smooth surface in M, called conformal boundary 
J. Physically it represent "infinity". In [15J, it is shown, that conformal boundaries 
must be spacelike hypersurfaces with respect to the conformal metric for all solutions of 
Eq. (fl4"|) with A > 0. One calls such solutions "future asymptotically de-Sitter" (FAdS) 
|13[ H], if its conformal boundary has a smooth non-empty future component J + ; there 
is the analogous concept for the past time direction. In particular, the de-Sitter solution 
[22] is FAdS. Under these conditions, J + represents the infinite timelike future of M. 
Some of the asymptotic geometric properties of FAdS solutions are discussed in [3]. 

Friedrich introduced his conformal field equations (CFE), as reviewed for instance in 
|15j . in order to deal with conformally compactified solutions of Einstein's field equations. 
In these conformally invariant equations, the fundamental unknown is the conformal 
metric g and the conformal factor Vt related to the physical metric g. The non-trivial 
property of these equations is, that they are, first, equivalent to Einstein's field equations 
wherever Q > 0, and, second, yield regular hyperbolic evolution equations even where 
Q = 0. Under the assumptions above, the CFE allow us to formulate what we call 
Cauchy problem" [13] . The idea is to prescribe data for the CFE on the hypersurface 
including its manifold structure, subject to certain constraints implied by the CFE. 
These data can then be integrated into the past by means of evolution equations implied 



16 



by the CFE. Friedrich proved that the J^-Cauchy problem is well-posed, and that the 
unique FAdS solution corresponding to a given choice of smooth data on J + is smooth, 
as long as it can be extended into the past. It is remarkable that this Cauchy problem 
allows to control the future asymptotics of the solutions explicitly by the choice of the 
data on J + . Concerning the past behavior of the solution corresponding to a given 
choice of data on however, there is only limited understanding and a-priori control, 
because of the complexity of the field equations. In this paper, we will give no details on 
the constraints on and say only briefly what the relevant initial data components 
are, since we do not want to introduce all necessary geometric concepts now. However, 
we write down a special class of solutions of the constraints in Section 13.1.21 We refer 
to [HIS], where the details have been carried out. 

We decided to use the so-called general conformal field equations, which are the 
CFE in conformal Gauss gauge [HI H5]. In our applications, we specialize the gauge 
even further to what we call Levi-Civita conformal Gauss gauge [3]. In this paper here, 
we will discuss neither the physical properties, nor the possibly bad implications of this 
choice of gauge, but just refer to [21 Hj. In any case, assuming, without loss of generality, 
A = 3, and having fixed the residual gauge initial data, as described in [3], the implied 
set of evolution equations is 





dte a c = 


Xa &b ' 


(15a) 




dtXab = 


-XaXcb - &E ab + L ab , 


(15b) 




dtT a \ = 


-Xa T dc + ^B ad e b c d , 


(15c) 




dtL a b = 


—d t U E ab — x a c L c b, 


(15d) 


d t E fe 


- D e c B a{f e ac e) = 


-ZXcEfe + 3X(e E f)c ~ XcEb C 9ef, 


(15e) 


d t B fe 


+ D ec E a(f£ aC e ) = 


-2>Xc C Bfe + 3X( e c -B/) c - X c b Bb C 9ef, 


(15f) 




fi(t) = 




(15g) 



for the unknowns 

u = ( e a 6 ; Xab, F a b c , L ab , E fe , B fe ^j . (15h) 

The unknowns are the spatial components e a b of a smooth frame field {e^} as in Eq. (jllbj) . 
with eo = dt, which is orthonormal with respect to the conformal metric, the spatial 
frame components of the second fundamental form \ab of the t = const-hypersurfaces 
with respect to eo, the spatial connection coefficients r a b c , given by T a b c e b = V ea e c — Xac&o 
where V is the Levi-Civita covariant derivative operator of the conformal metric, the 
spatial frame components of the Schouton tensor L ab , which is related to the Ricci tensor 
of the conformal metric by 

Lij = Rij/2 — gijg kl Rki/12, 

and the spatial frame components of the electric and magnetic parts of the rescaled 
conformal Weyl tensor E ab and B ab [T5J [16], defined with respect to eo- In this special 
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conformal Gauss gauge, the timelike frame field eo is hypersurface orthogonal, i.e. (Xab) 
is a symmetric matrix. In order to avoid confusions, we point out that, in principle, 
the conformal factor f2 is part of the unknowns in Friedrich's formulation of the CFE. 
However, for vacuum with arbitrary A, it is possible to integrate its evolution equation in 
any conformal Gauss gauge explicitly [14|, so that O takes the explicit form Eq. Q15g) in 



our gauge. We note, furthermore, that, since (E a b) and (B a b) are tracefree by definition, 
we can get rid of one of the components for each of the two. Our simple minded choice is 
the 33-component by £33 = —En — E22', the same for the magnetic part. The evolution 
equations Eqs. (|15e|) and ()15f|) of E a b and B a b are derived from the Bianchi system [15]. 
In our gauge, the constraint equations implied by the Bianchi system take the form 

D ec E\ - e ab e B daXb d = 0, D ec B c e + e ab e E daXb d = 0. (16) 

Here, e abc is the totally antisymmetric symbol with e 12 3 = 1, and indices are shifted 
by means of the conformal metric. The other constraints of the full system above are 
equally important, but are ignored for the presentation here. Further discussions of that 
evolution system and the quantities involved can be found in the references above. Note 
that in Eqs. t|15ej> . (115fj) and (fTHl) . the fields {e a } are henceforth considered as differential 
operators, using Eq. (lllbl) and writing the fields {Y a } as differential operators in the 
Euler coordinate basis as in Eqs. ©. Seen as a system of partial differential equations, 
the system Eqs. (I15p is symmetric hyperbolic, and hence the initial value problem is 
well-posed. 

Note that in this gauge, our initial hypersurface J + corresponds to t = 0. The past 
conformal boundary, if it exists, corresponds to t = 2. Hence, our time coordinate runs 
backwards with respect to physical time. 

These equations hold without any symmetry assumptions. In the following we will 
assume that all unknown fields involved are Gowdy symmetric. For the 1 + 1-code, we 
need to derive the evolution equations of the matrix (T a a ) . The fact that the conformal 
metric g is 13-invariant, implies, after straight forward computations, that the matrix 
(T a b) '■= (T a a g a 'b) is antisymmetric. Our choice of frame transport is parallel transport 
with respect to the conformal metric. This implies, after some algebra, that 

dt(T a a ') = 0. 

Hence, since the original system of equations is symmetric hyperbolic, also the "indirect 
reduction" to 1 + 1-dimensions is symmetric hyperbolic. So, the initial value problem 
for these equations is also well-posed. 



3.1.2 A class of initial data 

As initial data on J + , we use the "Berger data", which are solutions of the constraints 
derived for U(l)- and Gowdy symmetry in [3] . Those data are close to data of the A-Taub- 
NUT solutions and hence are particularly interesting for the strong cosmic censorship 
conjecture [3]. Here, we restrict to Gowdy symmetry. Under the conventions above, 
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these data take the form 



(e a b ) =diag(l,l,a 3 ), (17a) 

(Xo6) = diag(-l,-l,-l), (17b) 

^2 = 0, r 1 2 3 = -l/a 3 , ^ = 0, ^3 = 1/03, (17c) 

r 2 2 3 = o, r 3 1 2 = i/a 3 -2a 3 , r 3 1 3 = o, r 3 2 3 = o, (i7d) 

{Lab) = diag((5 - 3/a 2 )/2, (5 - 3/a|)/2, (-3 + 5/a 2 )/2) , (17e) 

{B ab ) = diag(-4(l - 4)/al -4(1 - a|)/a|, 8(1 - a|)/a|), (17f) 

(E0 + C2W20 —y/2az C 2 Rew 2 i \ 

E + C 2 w 20 -V2a 3 C 2 lmw 2 i J. (17g) 

-\/2a 3 C 2 Rew2i -\/2a 3 C 2 Imu7 2 i -2(£ + C2 w 20 ) / 



The induced conformal 3- metric of J + is a Berger sphere with a free parameter a 3 > 0. 
The only inhomogeneous, i.e. space dependent part of the initial data is given by the 
components E a b- For our definition of the functions w np , consult [3J; we just note that, 
with respect to the Euler coordinates, we have 

W20 = cos x, W21 = sin x e~ tpl j\Fl. 

For all these data, one finds 

/ 2 0\ 

(T/) = -2 . 

\ 00/ 

In total this family of solutions of the constraints has three free parameters a 3 > 0, i?o € 
M and C 2 6 K. We remark that the reason for the strange names of these parameters is 
consistency with our notation in [3j. 

3.1.3 Boundary control for the 1 + 1-code 

In Section 12.4.21 we have motivated our boundary control approach for the 1 + 1-code. 
Because the analysis depends strongly on the particular equations and choice of frame 
transport, it was not possible to give a further discussion there in full generality. Hence, 
let us continue here for our special choice of equations and frame transport. Due to what 
was said before, we have 

/ 2 0\ 

(T/) = -2 , 

\ 00/ 

for all t, for our choice of initial data. In this case, we say, that the frame is "boundary 
adapted" . Now, we introduce the new fields 

E* n := {E xl + E 22 )/2, E* 22 := {E xl - E 22 )/2, (18) 
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and similar for the magnetic part B a b, so that the boundary system, introduced in 
Section 12.4.21 yields the following conditions at x = and X = 7r, 

E12 = Ei 3 = E22 = E23 = 0, B12 = Bi 3 = B22 = B23 = 0, (19) 

whereas E^ and B*^ are free. For all other symmetric invariant 2-tensor fields, for 
instance the 2nd fundamental form, we get analogous relations, but, in addition, the 3- 
3-components are free, if the tensor is not tracefree. Since the behavior of the connection 
coefficients T a b c can be derived on the symmetry axes as well, which are the only non- 
tensorial objects in our set of unknowns, we obtain a complete set of boundary conditions 
for all the unknowns. The quantity Norm^ 1 ^ is now defined as the sum of the actual 
absolute numerical boundary values of all those unknowns, which are supposed to be 
zero there according to these results. Monitoring Norm^^ in a numerical computation, 
yields information on how well the boundary conditions are satisfied. 

In order to implement the 1 + 1-code numerically, we write the unknowns in terms of 
the new electric and magnetic fields defined in Eq. (|18p . Actually, it would be better to 
introduce the analogous combinations of fields for the other unknowns, but this has not 
yet been done. Hence, so far, the code lacks a clean way of enforcing e.g. the boundary 
condition xu ~ X22 = at x = and ir. To circumvent this problem temporarily, we 
have decided to work with a "partial enforcement" scheme, which, at a given time of the 
evolution, enforces all boundary conditions except for those of this type. In addition, we 
monitor the quantity Norm^^, and so far this treatment has turned out to be sufficient. 

3.2 Numerical results for the test application 

In [3] , we have performed a couple of tests with the 2 + 1-code, discussed the findings 
and drew conclusions about the numerical method. Here, we rather focus on the 1 + 1- 
code, and show so far unpublished tests and discussions in Section 13.2. 11 Afterwards in 
Section 13.2.21 we also compare a simulation done with the 2 + 1-code and the 1 + 1-code 
directly. 

We just note that we have not made systematic investigations of the CFL condition 
for our codes yet. 

3.2.1 Analysis of computations with the 1 + 1-code 

For our numerical test case here, we choose 03 = 0.7, C2 = 0.1 and Eq = in Eqs. (|17p 
as initial data parameters, corresponding to the "large inhomogeneity case" in |3] and to 
one of the simulations presented in [JJ. The associated solution turns out to develop a 
singularity, and hence can be seen as an interesting test case for our code. The evolution 
of a spatial norm of the curvature invariant called Kretschmann scalar is shown in Fig. [U 
All the results we show here were done without the automatic spatial adaption approach 
described in Section l2.4.1l because, in order to study convergence, it seems more useful to 
control and adapt the spatial resolution manually. The adaption norm, computed with 
respect to £13, was used only for estimating the aliasing error. The time integration 
was done with the adaptive 5th order embedded RK scheme with control parameters rj 
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and fo m in as discussed in Section 12.4.11 For these runs, we decided to use the "partial 
enforcement" scheme of the boundary conditions, explained in Section 13.1.31 All runs 
were done with double precision. 

The constraints Eq. (|16p are satisfied initially up to machine precision. However, 
due to numerical errors, those constraints typically become violated more and more 
with increasing evolution time. Let us define Norm^ cons * r ^ as the Z^-norm of the sum of 
the absolute values of each of the six components of the left hand sides of Eqs. (fT6|) at 
a given instant of time, all that divided by tr {xab)-> i n order to factor out the observed 
collapse of the solution. L p -norms of functions on § 3 are always evaluated here by means 
of their corresponding functions on T 3 and the standard L p -norm on T 3 . Another norm, 
which measures how well the numerical solution satisfies Einstein's field equations, is 



^ ovm ( einstein ) ■- 



L!(S 3 ) 



where the Ricci tensor Rij of the physical metric is evaluated algebraically from the 
conformal Schouton tensor Lij and derivatives of the conformal factor f2. The indices 
in this expression are defined with respect to the physical orthonormal frame given by 
Si = Qet. The norm is computed by summing over the Zi -norms of each component. 

Now, we will distinguish two phases of the evolution for these initial data, in which 
different aspects and effects are important: the early evolution close to J + for t between 
0.0 and 0.69, and the late evolution close to the singularity, which we find at t ~ 
0.69520493. In order to avoid confusions, we recall that the terms "early" and "late" are 
always understood with respect to the time coordinate t, which, however, runs backwards 
with respect to the physical time. 

At early times, it is achieved easily that the spatial discretization error is not signifi- 
cant, until some later time when small spatial structure starts to form more rapidly. One 
hint that this is true, as we do not show here, is that NomS adapi ^ is more or less constant 
over a long time period, and small. Another hint is that the behaviors of Norm^ consir - ) 
and Norm^ emstem - ) are not strongly influenced by the spatial resolution. See Fig. [2] and 
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Fig. 3: Violations of EFE at early times 
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Fig. 4: Violation of boundary conditions at 
early times 



[3l where N represents the number of spatial grid points, which is constant in this early 
regime. Indeed, the higher N, the larger is the initial value Norm^ cons * r ^ due to higher 
round-off errors for computing spatial derivatives. This is not visible for Norm^ em,s * em \ 
since this quantity is defined purely algebraically in the unknowns. In Fig. EJ we see 
that Norm( co ™* r ) grows less, the higher the time resolution is, i.e., in particular, the 
smaller the parameters rj and eventually also h m i n are. However, we always observe 
at least weak approximately exponential growth. In Fig. El we see a similar behavior 
for Norm( ems ' em ). We do not show here that there is neither a particular growth of 
Norm( cons4r ) , nor of Norm^ emstem \ at the symmetry axes. Rather, the maximal growth 
takes place, where the curvature increases most strongly. This can be seen as a confir- 
mation that our treatment of the coordinate singularities works well, cf. [3]. Note that 
there is an optimal time resolution, in the sense that, if we choose a higher resolution, 
the constraint error and Norm^ ems * em ^ are actually increased caused by higher round-off 
errors. Although we do not show any plots, we want to mention, that we have experi- 
mented with "quad precision". For the 1 + 1-code, this yields reasonable performance 
and has several consequences. First, the initial data for the constraint violations are 
decreased by many orders of magnitude, since those are determined primarily by the 
precision of the numerical number representation. By choosing appropriate resolutions, 
we find that the constraint violations and Norm^ ems * em ^ can then be kept several orders 
of magnitude smaller than in the double precision case during the whole run. However, 
they always show exponential growth, which suggests that this is the typical behavior 
of the constraint propagation in our system of equations. Furthermore, quad precision 
allows us to work in a regime in which discretization errors are much larger than round- 
off errors, and hence it is easier to interpret convergence tests. Fig. 0] shows the behavior 
ofNorm^. The errors at the boundaries are small and stable, despite of some weak 
growth, which is expected in situations close to a singularity. In these tests, this error 
does not decrease for higher resolutions, which is a hint that round-off errors have a 
significant effect. With quad precision, we were able to confirm that the violations of 
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Fig. 5: Adaption norm at late times Fig. 6: Constraint violations at late times 



the boundary conditions become smaller, consistent with increasing resolution. 

Concerning the late evolution, the following is found. The following figures, Fig. [5] 
to Fig. show a very small time neighborhood of the final time, where the runs were 
stopped and where the solution blows up. In the runs underlying these plots, we adapted 
the spatial resolution several times during the runs manually; the number N in the figures 
is the final spatial resolution in each case. Each manual spatial adaption step is visible 
in the plots as a jump, because for different spatial resolutions, the numerical values of 
the norms slightly change. In all these runs, we choose rj = 10 -13 , and hence the time 
steps h decrease so strongly that h = h m i n at that time when the runs were stopped. In 
this late time regime, the errors are dominated by spatial discretization errors, because 
the solution has the property that spatial structures shrink without bound. This can 
be seen by looking at the late time plot of the adaption norm in Fig. O It shows, how 
strongly the demand for spatial resolution grows with time, but also, that it is possible to 
gain control by increasing the resolution at least temporarily. However, the demand for 
spatial resolution increases very strongly with time, and it turns into a difficult numerical 
issue to keep track of that eventually. In Fig. El we demonstrate, how the choice of spatial 
resolution influences the propagation of the constraint violations, and that this quantity 
converges to a weakly exponentially growing for sufficiently high spatial resolutions. This 
is a promising result, and shows that the constraint propagation is more or less under 
control, as long as it is possible to increase the resolution in practice. Fig. \7\ shows the 
violations of the boundary conditions. They turn out to be very small and under control 
for sufficiently large resolutions. The higher the final spatial resolution is, the smaller 
these violations appear. Finally, in Fig. [8] we show Norm( ems * em ). We do not observe 
a very strong difference between the various resolutions; indeed, in order to make the 
differences visible at all, the time axis represents an even smaller time neighborhood 
now. It is unexpected that at very late times this norm is not necessarily smaller the 
higher the resolution is, and it has to be investigated whether this is a problem. 
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3.2.2 Comparison of a computation with the 2 + 1- and 1 + 1-code for a 
Gowdy symmetric solution 

Next, we want to compare directly the numerical results obtained with the 2 + 1-code 
and the 1 + 1-code for a Gowdy symmetric solution of our equations. We restrict to a 
singularity free case, since we do not want to be spoiled by the lack of spatial resolution, in 
particular for the 2+ 1-code. The following initial data parameters are chosen: 03 = 0.93, 
Eq = and C2 = 0.5, which correspond to the "regular case" in [3]. The corresponding 
solution is smooth for all < t < 2 and hence develops a smooth past conformal 
boundary at t = 2. Fig. shows the evolution of the curvature invariant Kretschmann 
scalar. We show neither convergence plots, nor the error quantities before, since the 
situation is very similar to the early phase of the singular solution in the previous section. 
Instead, we report only on one run done with one fixed resolution: the size of the time 
step is h = 5 • 10~ 4 and the number of spatial points is N = 40 for the 1 + 1-code and 
N x = 41, N pi = 21 for the 2 + 1-code. We use the non-adaptive 4th order RK time 
integrator, and the automatic spatial adaption has been switched off here as well. 
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Fig. 10: Violations of the boundary condi- 
tions of the 1 + 1-code 
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Fig. 11: Comparison of the 2 + 1- and 1 + 1- 
code 



The run with the 1 + 1-code was done without enforcing the boundary conditions 
now, cf. Fig. [TUl We see that this error stays very small and is stable. In order to 
illustrate, how well the results obtained with the 1 + 1- and the 2 + 1-code coincide, let 
us define the following norm 



Norm^ :-- 



(i+i) I 



\X =7T J 



(2+1) I 



| x =7r,pi=0 y 



Consider Fig. [TT] where this norm is plotted vs. time. We find very good agreement 
between the two codes. Some deviations can be expected, since in the 2 + 1-code, Gowdy 
symmetry is valid only approximately; see more comments in [3j. So far, we have made 
no effort to explain the oscillatory behavior in both figures, but we conjecture them to be 
caused by aliasing, since its amplitude becomes smaller, the higher the spatial resolution 
is. We have also compared more variables at other grid points and found similar results. 



4 Summary and conclusions 

The purpose of this paper is to introduce our numerical approach for the study of 
evolution equations with spatial §> 3 -topology. First, we have given details on a geometric 
point of view in order to find a natural discretization. Second, we have discussed issues 
related to the implementation and, third, analyzed test applications. By all this we 
have demonstrated the feasibility of this approach. Indeed, these techniques have been 
applied elsewhere, for instance in [21 El [5]. 

Although our main interest lies in the field of general relativity, we believe that the 
applicability of the method is more general. For the presentation of this paper, we made 
a couple of special choices which need to be overcome in order to use the method in more 
general situations. First, our analysis is based on U(l)-symmetry so far. We believe that 
it is straight forward to generalize in principle. Furthermore, we decided to formulate 
the equations in terms of smooth global frames on S 3 . This had the consequence that 
only singular terms of special type are present, and with this knowledge, we were able 
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to work out the Fourier series of the formally singular terms explicitly in Section 12, 21 
However, in general, as soon as we have a-priori knowledge about the structure of the 
formally singular terms in the equations, and if their type is not completely different 
than ours, a similar regularization is possible. Hence, from this point of view, one should 
be able to modify our method to cases, where other kinds of frames or spatial coordinates 
are used. In our application, we made a special choice of coordinate gauge and frame 
transport. This implied that, in particular, the matrix (T a a ) is constant in time. If the 
original evolution equations imply a well-posed initial value problem, then, in this gauge, 
also the evolution equations for Gowdy symmetry reduced to 1 + 1 dimensions have this 
property. For other gauges, however, this might not be the case. A further issue turns 
up when the frame transport does not keep the orthonormal frame boundary adapted in 
the sense of Section \3. 1.31 Then a different boundary treatment than ours might become 
necessary. As a last point of this certainly incomplete list, we mention that we have 
not made experiments with other time integrators yet. Of particular interest might be 
implicit schemes, for instance, to treat parabolic evolution equations. 

We believe that our current numerical infrastructure has not yet been pushed to its 
limits, but is also not yet really optimized for high spatial resolutions. For instance, 
we still do not use FFT but only the partial summation scheme. Also it may be true, 
that there is a more optimal trade-off between accuracy and efficiency for other time 
integrators than the Runge Kutta schemes of our choice; comments on this can be found 
in [Hj. Furthermore, it might make sense to think about parallelization of the code. This 
should be straight forward with some publicly available FFT libraries, for instance |17j . 
We expect, that the 1 + 1-code, in its current form, is limited to a few thousand spatial 
grid points. We want to stress that this turned out to be sufficient in our runs with T 3 - 
topology in [3]. There, we were able to reproduce "spiky features" numerically. Hence, 
we are optimistic that also such difficult phenomena can be studied with our numerical 
infrastructure in 1 + 1 dimensions. Dependent on the kind of the application, one may 
nevertheless doubt, if spectral methods are suitable at all. In general, it is fair to say 
that, although these methods are highly accurate for lower resolutions, they might be 
too inefficient for high resolutions. Thus it makes sense, to also investigate into other 
methods, for instance finite differencing methods, maybe with multipatch [38| . II 1 |, [27] . in 
order to make systematic comparison studies. Further automatic local spatial adaption, 
usually called adaptive mesh refinement, would be desirable; for instance cosmological 
solutions with spikes have been investigated with such techniques in [20J, however not 
for § 3 -topology. 

In summary, due to the results of the tests here and in O [U [5], we believe that 
our numerical approach has promising potentials, which have not been fully exhausted 
yet, for many kinds of applications with spatial § 3 -topology and beyond. For instance, 
we find it particularly appealing that our infrastructure can be used directly also for 
problems with spatial T 3 -topology and, at least for Gowdy symmetry, with S 1 x S 2 - 
topology [5]. For the applications, which we have in mind, future research will show, 
how much the method must be modified, or whether completely new approaches will 
become necessary, in order to deal with probably much more severe phenomena than 
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those present in the classes of solutions considered so far. 
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